Pré-requisitos

Pacotes necessários:

#if(!require(PSF)) install.packages("PSF")
#if(!require(timetk)) remotes::install_github("business-science/timetk")
if(!require(EflowStats)) remotes::install_github("USGS-R/EflowStats")


pacotes <- c(
  "here",
  "usethis",
  "data.table",
  "HEobs",
#  "PSF",
  "tidyverse",
  "lubridate",
  "fs",
  "checkmate",
#  "xts",
#  "hydroGOF",
#  "ModelMetrics",
#  "forecast",
#  "timetk",
  "EflowStats",
  "hydrostats",
 #"NbClust",
 #  "cluster",  
 #  "cluster.datasets", 
  "cowplot", 
  # "clValid",
  "ggfortify", 
  #"clustree",
  #"dendextend",
  #"factoextra",
  #"FactoMineR",
  #"corrplot",
  #"GGally",
  #"ggiraphExtra",
  "kableExtra",
  "tidymodels",
  "fasstr"
)
# Carregar os pacotes
easypackages::libraries(pacotes)

Scripts:

source(here('R', 'load-data.R'))
source(here('R', 'utils.R'))
source(here('R', 'prep-stn-data.R'))

Dados de vazão

Os dados importados de vazão devem ser regularmente espaçados no tempo. Esta adequação das séries diárias, se necessária, pode ser realizada com a função complete_dates() do pacote {lhmetools}. Assim assegura-se que todas estações possuam 1 observação por dia e sem datas faltantes.

Metadados

data_link <- "https://www.dropbox.com/s/d40adhw66uwueet/VazoesNaturaisONS_D_87UHEsDirceuAssis_2018.dat?dl=1"
qnat_meta <- extract_metadata(NA_character_, informative = TRUE)
glimpse(qnat_meta)
## Rows: 87
## Columns: 5
## $ estacao_codigo <dbl> 18, 237, 215, 119, 190, 32, 14, 247, 1, 216, 191, 149, …
## $ latitude       <dbl> -19.918751, -22.636078, -27.776389, -23.811460, -6.7985…
## $ longitude      <dbl> -49.92053, -48.27274, -51.18944, -46.52390, -43.92756, …
## $ nome_estacao   <chr> "A. VERMELHA", "BARRA BONITA", "BARRA GRANDE", "BILLING…
## $ municipio      <chr> "A. VERMELHA", "BARRA BONITA", "BARRA GRANDE", "BILLING…

Dados

qnat_data <- qnat_dly_ons() %>%
  select(date, qnat, code_stn) %>%
  lhmetools::complete_dates(group = "code_stn")
glimpse(qnat_data)
## Rows: 2,796,267
## Columns: 3
## $ date     <date> 1931-01-02, 1931-01-03, 1931-01-04, 1931-01-05, 1931-01-06, …
## $ code_stn <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ qnat     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
# Incluindo nome da estacao
qnat_data <- qnat_data %>%
  full_join(select(qnat_meta, estacao_codigo, nome_estacao),
            by = c(code_stn = "estacao_codigo")
            ) %>%
  rename("name_stn" = nome_estacao)

Assinaturas hidrológicas para 1 posto

Preparação dos dados para aplicação da função calc_magnifSeven com a opção de ano hidrológico (water year) para o posto 74 da ONS (G. B. Munhoz).

qnat_posto <- qnat_data %>% 
  sel_station(.,station = 74)  
glimpse(qnat_posto)
## Rows: 32,141
## Columns: 4
## $ date     <date> 1931-01-02, 1931-01-03, 1931-01-04, 1931-01-05, 1931-01-06, …
## $ code_stn <dbl> 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 74, 7…
## $ qnat     <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
## $ name_stn <chr> "G.B. MUNHOZ", "G.B. MUNHOZ", "G.B. MUNHOZ", "G.B. MUNHOZ", "…
summary(qnat_posto)
##       date               code_stn       qnat           name_stn        
##  Min.   :1931-01-02   Min.   :74   Min.   :  69.58   Length:32141      
##  1st Qu.:1953-01-01   1st Qu.:74   1st Qu.: 300.02   Class :character  
##  Median :1975-01-01   Median :74   Median : 524.20   Mode  :character  
##  Mean   :1975-01-01   Mean   :74   Mean   : 730.61                     
##  3rd Qu.:1996-12-31   3rd Qu.:74   3rd Qu.: 958.57                     
##  Max.   :2018-12-31   Max.   :74   Max.   :7832.77                     
##                                    NA's   :13605
magnif7(stn_data = select(qnat_posto, date, qnat))
## # A tibble: 8 × 2
##   indice    statistic
##   <chr>         <dbl>
## 1 lam1        736.   
## 2 tau2          0.42 
## 3 tau3          0.35 
## 4 tau4          0.18 
## 5 ar1           0.982
## 6 amplitude     0.17 
## 7 phase       262.   
## 8 bfi           0.413

Dados agrupados por postos

by_stn <- qnat_data %>% 
  group_by(code_stn) %>%
  nest()

Testando fasstr

https://bcgov.github.io/fasstr/articles/fasstr.html

Será relevante considerar a diferenças nos períodos dos anos hidrológicos por posto?

check_season_plots <- qnat_data %>%
  rename("Date" = date) %>%
  plot_daily_stats(
    values = "qnat",
    groups = "code_stn",
    start_year = 1991,
    end_year = 2010,
    log_discharge = TRUE,
    add_year = 2001, 
    complete_years = TRUE, 
    include_title = TRUE
  )
codes <- names(check_season_plots) %>%
  readr::parse_number()

lookup <- qnat_data %>% 
  select(contains("stn")) %>% 
  distinct() %>%
  slice(order(codes)) %>%
  arrange(code_stn)

plot_l <- map(seq_along(codes), 
    function(i) {
      # i = 1
      nm <- filter(lookup, code_stn == codes[i]) %>%
        pull(name_stn)
      check_season_plots[[i]] + ggtitle(paste0("Posto: ", nm))
    })
plot_l
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

## 
## [[18]]

## 
## [[19]]

## 
## [[20]]

## 
## [[21]]

## 
## [[22]]

## 
## [[23]]

## 
## [[24]]

## 
## [[25]]

## 
## [[26]]

## 
## [[27]]

## 
## [[28]]

## 
## [[29]]

## 
## [[30]]

## 
## [[31]]

## 
## [[32]]

## 
## [[33]]

## 
## [[34]]

## 
## [[35]]

## 
## [[36]]

## 
## [[37]]

## 
## [[38]]

## 
## [[39]]

## 
## [[40]]

## 
## [[41]]

## 
## [[42]]

## 
## [[43]]

## 
## [[44]]

## 
## [[45]]

## 
## [[46]]
## Warning: Transformation introduced infinite values in continuous y-axis

## 
## [[47]]

## 
## [[48]]

## 
## [[49]]

## 
## [[50]]

## 
## [[51]]

## 
## [[52]]

## 
## [[53]]

## 
## [[54]]

## 
## [[55]]

## 
## [[56]]

## 
## [[57]]

## 
## [[58]]

## 
## [[59]]

## 
## [[60]]

## 
## [[61]]

## 
## [[62]]

## 
## [[63]]

## 
## [[64]]

## 
## [[65]]

## 
## [[66]]

## 
## [[67]]

## 
## [[68]]

## 
## [[69]]

## 
## [[70]]

## 
## [[71]]

## 
## [[72]]

## 
## [[73]]

## 
## [[74]]

## 
## [[75]]

## 
## [[76]]

## 
## [[77]]

## 
## [[78]]
## Warning: Transformation introduced infinite values in continuous y-axis

## 
## [[79]]

## 
## [[80]]

## 
## [[81]]

## 
## [[82]]

## 
## [[83]]

## 
## [[84]]

## 
## [[85]]

## 
## [[86]]

## 
## [[87]]

Assinaturas hidrológicas para todos postos

seven_stats <- by_stn %>%
  ungroup() %>%
  mutate(stats = map(data, ~.x %>% select(date, qnat) %>%magnif7))
seven_stats_tidy <- seven_stats %>%
  select(stats) %>%
  unnest() %>%
  pivot_wider(names_from = indice, values_from = statistic) %>%
  ungroup()
## Warning: `cols` is now required when using unnest().
## Please use `cols = c(stats)`
## Warning: Values are not uniquely identified; output will contain list-cols.
## * Use `values_fn = list` to suppress this warning.
## * Use `values_fn = length` to identify where the duplicates arise
## * Use `values_fn = {summary_fun}` to summarise duplicates
seven_stats_tidy
## # A tibble: 1 × 8
##   lam1       tau2       tau3       tau4       ar1        amplitude  phase  bfi  
##   <list>     <list>     <list>     <list>     <list>     <list>     <list> <lis>
## 1 <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl [87]> <dbl … <dbl…
saveRDS(seven_stats_tidy, file = here("output", "seven_stats.RDS"))